'''
Created on Aug 20, 2009

@author: echeek
'''
import math
import primesieve

def Solovay(n):
    import random
    for k in range(10):
         a=random.randint(1,n-1);
         x=a/n;
         if x==0 or (a**((n-1)/2) % n)!=x:
             return False;
    return True;

def Euler216():
    maxroot=50000000#70710679;
    '''
    sieve=primesieve.bitSieve(maxroot);
    print 'sieve generated';
    i=3;
    mark=3;
    primes=[2];
    while i<maxroot:
       if (mark & sieve)==0:
           primes.append(i);
           #print i;
       mark<<=2;
       i+=2; 
    print 'primes copied';
    '''
    primes=primesieve.listPrimesBasic(maxroot)
    print 'primes list complete';
    count=0;
    n=2;
    while n<=50000000:
        tn=2*(n**2)-1;
        stoppoint=math.sqrt(tn);
        isprime=Solovay(tn);
        if(isprime):
            for p in primes:
                if(p>stoppoint):
                    break;
                if( tn % p==0):
                    isprime=False;
                    break;
        if(isprime):
            #print "%i:%i" % (n,tn);
            count+=1;
        n+=1;
        if n%100==0:
            print n;
    print count;
    

             
    
def Euler216b():
    maxroot=70710679;
    '''
    sieve=primesieve.bitSieve(maxroot);
    print 'sieve generated';
    i=3;
    mark=3;
    primes=[2];
    while i<maxroot:
       if (mark & sieve)==0:
           primes.append(i);
           #print i;
       mark<<=2;
       i+=2; 
    print 'primes copied';
    '''
    primes=primesieve.listPrimesBasic(maxroot)
    
    
Euler216();
